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We utilize a) detailed end-to-end numerical simulations of sample narrow-angle as- 
' trometric observing campaigns with the Space Interferometry Mission (SIM) and the 

subsequent data analysis process, and b) the set of extra-solar planetary systems dis- 
covered so far by radial velocity surveys as templates to provide meaningful estimates 
of the limiting capabilities of SIM for the detection and measurement of multiple-planet 
systems around solar-type stars in the solar neighborhood. 
.£h ! We employ standard \ 2 statistics, periodogram and Fourier analysis to evaluate 

SIM's ability to detect multiple planetary signatures: the probability of detecting ad- 
ditional companions is essentially unchanged with respect to the single-planet config- 
urations, but after fitting and subtraction of orbits with astrometric signal-to-noise 
ratio a/ ad — * 1 the false detection rates can be enhanced by up to a factor 2; the 
periodogram approach results in robust multiple-planet detection for systems with pe- 
riods shorter than the SIM mission length, even at low values of a/ ad, while the least 
squares technique combined with Fourier series expansions is arguably preferable in the 
long-period regime. We explore the three-dimensional parameter space defined by as- 
trometric signature, orbital period, and eccentricity to derive general conclusions on the 
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capability of SIM to accurately measure the full set of orbital parameters and masses 
for a variety of configurations of planetary systems: the accuracy on multiple-planet 
orbit reconstruction and mass determination suffers a typical degradation of 30-40% 
with respect to single-planet solutions; mass and orbital inclination can be measured 
to better than 10% for periods as short as 0.1 yr, and for a/ad as low as ~ 5, while 
ot/ud — 100 is required in order to measure with similar accuracy systems harboring 
objects with periods as long as three times the mission duration. We gauge the potential 
of SIM for meaningful coplanarity measurements via determination of the true geome- 
try of multiple-planet orbits: for systems with all components producing a /ad — 10 or 
greater, quasi-coplanarity can be reliably established with uncertainties of a few degrees, 
for periods in the range 0.1 < T < 15 yr; in systems where at least one component has 
a / CJ d — * 1) coplanarity measurements are compromised, with typical uncertainties on 
the mutual inclinations of order of 30° — 40° . We quantify the improvement derived in 
full-orbit reconstruction and planet mass determination by constraining the multiple- 
planet orbital fits to SIM observations with the nominal orbital elements obtained from 
the radial velocity measurements: the uncertainties on orbital elements and masses can 
be reduced by up to an order of magnitude, especially for long-period orbits in face-on 
configurations, and for low amplitude orbits seen edge-on. 

Our findings are illustrative of the importance of the contribution SIM will make, 
complementing other on-going and planned spectroscopic, astrometric, and photometric 
surveys, in order to fulfill the expectations for ground-breaking science in the fields of 
formation and evolution of planetary systems during the next decade. 

Subject headings: astrometry - planetary systems - instrumentation: interferometers - 
methods: data analysis - methods: numerical 

1. Introduction 

After more than thirteen years of precise radial velocity measurements, today there is clear 
evidence that several nearby solar-type stars harbor candidate planetary systems, composed of 
two or more planets (Butler et al. 1999; Marcy et al. 2001b; Mayor et al. 2001; Jones et al. 2002; 
Fischer et al. 2002; Marcy et al. 2002b; Fischer et al. 2003; Butler et al. 2003), or a planet and a 
probable brown dwarf (Marcy et al. 2001a; Udry et al. 2002) 

The answers provided by these discoveries have presented even more challenging puzzles to 
theoretical models aimed at describing the formation and evolution of planetary systems. The va- 
riety of orbital arrangements in multiple-planet systems found by spectroscopy calls into question 
both the origin and early dynamical evolution of such systems, in terms of formation mecha- 
nisms and orbital migration scenarios, as well as their long-term dynamical evolution and sta- 
bility. Early attempts have been made to connect and relate the early stages of their formation 
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to their presently observed orbital properties. For example, the two planets orbiting the M4V 
star Gliese 876 appear to be locked in a 2:1 resonance, and the existence of such commensura- 
bilities may be considered important evidence of possible inward migration for pairs of planets 
due to mutual interactions as well as tidal interaction with the protoplanetary disk during and 
shortly post-formation (Snellgrove et al. 2001; Nelson & Papaloizou 2002). Similarly, the three- 
planet system around the F8V star v Andromeda? (v And) appears to be stable over evolution 
timescales of order of the star's lifetime, again suggesting the possibility that the system's or- 
bits were shaped by the interaction of protoplanets with the protoplanetary disk, in the first few 
million years of v And's existence (Artymowicz 2001). Furthermore, in many of the observed 
multiple-planet systems the eccentricities are much larger than those of the planets of our own 
solar system. Again, this may be explained with significant eccentricity evolution due to mutual 
interaction between planets as well as with the protoplanetary disk during the epoch of disk de- 
pletion (Chiang & Murray 2002; Nagasawa et al. 2003). However, direct extrapolations from the 
early stages of their formation to the presently observed configurations of extra-solar planetary 
systems are still somewhat speculative. Today, in fact, there is still a lack of observational support 
to the theoretical models describing the early epochs of the formation of planetary systems, either 
because of the objective difficulty to identify good targets among very young stellar objects (often 
completely obscured by circumstellar material) or due to insufficient sensitivity of the present gener- 
ation of ground-based and space-borne instrumentation. Thus, theorists have focused primarily on 
the issues of long-term stability and orbital evolution of extra-solar planetary systems, independent 
on the details of their formation. 

Since Newton's discovery of the law of universal gravitation, the problem of predicting the 
orbits of the planets in our own solar system has been tackled by many giants of mathematics 
(Euler, Laplace, Jacobi, Lagrange, Gauss, Poincare). By the end of the nineteenth century, it had 
been demonstrated that it is not possible to find an exact analytic solution describing the motion 
of iV bodies, for A > 3, and the theory of perturbations in celestial mechanics was developed and 
extensively applied to many interesting cases, always within the boundaries of our solar system. 
During the two last decades, with the advent and exponential progress of computer technology, 
large-scale computations that follow the full evolution of the nine planets for a significant fraction 
of the lifetime of the Sun became possible. This has allowed the establishment of, for example, the 
long-term chaotic behavior of the orbits of planets in the solar system (Sussman & Wisdom 1992). 
In fact, the typical (Lyapunov) timescale for exponential divergence of two orbits in the solar system 
starting with slightly different initial conditions is only several million years, and soon the issue of 
stability only becomes a statistical question. The planetary orbits in the solar system are on the 
other hand stable, in the sense that the probability of ejections or close encounters during a time 
frame comparable to the main-sequence lifetime of the Sun (~ 12 Gyr) is extremely low. 

However, until the first multiple-planet system was announced (Butler et al. 1999), the prob- 
lem of planetary orbits dynamics had in recent years been almost relegated to the level of an 
academic pursuit. Now that a variety of additional systems, with a wide range of different orbital 
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arrangements, has been found, the field of planetary system dynamics has received revived inter- 
est. The long-term dynamical evolution of the 11 extra-solar planetary systems known to date has 
been extensively investigated by theorists, utilizing both direct numerical integrations and approx- 
imate analytical approaches for the description of the N-body mutual interactions. As opposed to 
our own solar system, or the three-planet system discovered by Wolszczan k, Frail (1992) around 
the pulsar PSR 1257+12, both appearing dynamically stable over long timescales (Laskar 1994; 
Gladman 1993), the dynamical evolution and long-term stability issues for the majority of the 
multiple-planet systems discovered by radial velocity are highly uncertain. As a general result, 
studies specifically targeted to gauge the general dynamical behavior of extra-solar planetary 
systems such as v And (Laughlin & Adams 1999; Rivera & Lissauer 2000; Stepinski et al. 2000; 
Jiang & Ip 2001; Barnes & Quinn 2001; Chiang et al. 2001; Rivera & Lissauer 2001a, Chiang & 
Murray 2002), 55 Cancri (55 Cnc hereafter. Novak et al. 2002; Ji et al. 2003), 47 Ursae Majoris (47 
UMa hereafter. Godzdziewski 2002; Laughlin et al. 2002), Gliese 876 (Kinoshita & Nakai 2001; 
Laughlin & Chambers 2001; Rivera & Lissauer 2001b; Ji et al. 2002a; Godzdziewski et al. 2002; 
Lee & Peale 2002a, 2002b), HD 82943, HD 37124, HD 12661, HD 38529, and HD 160691 (Godzdziewski 
& Maciejewski 2001; Kiseleva-Eggleton et al. 2002; Godzdziewski & Maciejewski 2003; Godzdziewski 
2003a, 2003b), come to the same conclusion: the stability of the systems can be greatly affected 
by small variations of the instantaneous orbital elements of the planetary orbits obtained from 
the radial velocity data and utilized as initial conditions for the numerical integrations. In par- 
ticular, for the systems not to be destabilized and disrupted on very short timescales (as short 
as 10 3 -10 5 years), constraints must be placed on the maximum allowed values for the masses and 
on the range of allowed relative inclinations of the orbits. Often analytical calculations are per- 
formed with the coplanarity of the orbits as a working assumption. Early attempts have also 
been made to verify the possibility of regions of dynamical stability inside the parent stars' Hab- 
itable Zones (Kasting et al. 1993), where Earth-size planets may be found (Jones et al. 2001; 
Noble et al. 2002; Jones & Sleep 2002a, 2002b; Rivera & Haghighipour 2002; Cuntz et al. 2003; 
Menou k, Tabachnik 2003). However, the existence of such regions does not directly imply that 
rocky planets may have actually formed, in these systems, at such privileged distances, in the first 
place (Thebault et al. 2002). 

The results obtained by many of the present studies on the dynamical stability of extra-solar 
planetary systems contain significant ambiguities. The impossibility to draw at present more than 
general statistical conclusions on such issues arises partly from the intrinsically chaotic nature of 
N-body systems, and partly from a lack of knowledge of a few key parameters that cannot be 
derived by the present observational datasets, but must be a priori fixed throughout the analyses. 
As a matter of fact, in general four quantities drive the effects of the mutual perturbations in 
multiple-planet systems investigated by theorists: the planetary eccentricities, the orbital periods 
(and possible commensurabilities among them), the planet masses (usually through their ratio), 
and the relative inclinations of the orbital planes. The last two of these parameters cannot be 
determined from the intrinsically one-dimensional information extracted from stellar spectra. In 
fact, radial velocity measurements do not determine either the inclination i of the orbital plane with 
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respect to the plane of the sky (and this in turn allows only for lower limits to be placed on the 
actual mass of each companion to the observed star) or the position angle f2 of the line of nodes in 
the plane of the sky (and thus relative inclination angles remain unknown): without knowledge of 
the full three-dimensional geometry of the systems and true mass values, general conclusions on the 
architecture, orbital evolution and long-term stability of the newly discovered planetary systems 
remain questionable. 

Attempts have been made to break the sini degeneracy and determine the true masses of 
planets in multiple systems utilizing self-consistent dynamical fits (Laughlin & Chambers 2001; 
Rivera & Lissauer 2001b). Such procedures significantly improve orbital fits performed assuming 
independent Keplerian motions in cases such as that of Gliese 876, where mutual perturbations 
are relevant and the orbital elements of the planets undergo significant changes on timescales 
comparable to the timespan of the observations. However, these dynamical fits to the radial velocity 
data are not conclusive, as assumptions must still be made on the actual relative inclinations of the 
planets, and broad ranges of sini provide similarly good fits. Thus, tightly constraining planetary 
masses will be difficult for several years to come, if one has to rely only upon radial velocity 
measurements. Furthermore, different orbital configurations giving similarly good values of the 
reduced chi-square xt = X 2 l v (where v is the number of degrees of freedom) behave very differently 
when integrated over long timescales, and the ambiguity on the actual long-term stability of the 
systems cannot be removed. 

In the latest years other techniques have reached the sensitivity necessary to complement 
spectroscopy in extra-solar planet searches. In particular, transit photometry has been success- 
ful in identifying the first gas giant planet eclipsing its parent star (Charbonneau et al. 2000; 
Henry et al. 2000), confirming the spectroscopic detection (Mazeh et al. 2000; Queloz et al. 2000), 
and very recently in directly discovering the second transiting planet (Konacki et al. 2003). The 
observations of planetary transits on the disk of HD 209458 allow for a direct estimate of the planet's 
size and actual mass (as the inclination of the orbital plane is known, and the sin i degeneracy breaks 
down), thus, also thanks to the recent detection of sodium in its atmosphere (Charbonneau et al. 2002), 
its density can be inferred and important insights on its composition can be obtained. How- 
ever, transit photometry detects only the small proportion of planets whose orbits happen to 
line up almost exactly with the line of sight to the star. We anticipate that future high-precision 
ground-based (Mariotti et al. 1998; Booth et al. 1999; Colavita et al. 1999) as well as space-borne 
(Danner & Unwin 1999; Perryman et al. 2001) astrometric observatories will be among the most 
effective techniques to remove the present ambiguities in the dynamical analysis of extra-solar 
planetary systems as well as provide valuable data to probe models addressing the issues of the 
formation of planetary systems and the actual nature of sub-stellar companions. The intrinsically 
two-dimensional astrometric data provide the means to directly measure the two missing parame- 
ters, i.e. the inclination angle and the line of nodes of each planetary orbit. Knowledge of the full 
viewing geometry of a system of planets allows then to derive meaningful estimates for the true 
masses of each orbiting object and the relative inclinations of the orbital planes. 
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In our previous work (Sozzetti et al. 2002, S02 hereafter), we addressed the issues of the de- 
tectability and measurability of single planets around single stars with the Space Interferometry 
Mission (SIM), with the instrument operated in narrow-angle astrometric mode. We expressed our 
results as a function of both SIM mission parameters and properties of the planet (mass, orbital 
characteristics). In the continuation of our studies we have utilized extensive simulations of SIM 
observations of the present sample of extra-solar multiple-planet systems in order to draw general 
conclusions on the ability of SIM to discover and measure systems of planets, as well as quantify the 
instrument's capability to determine the coplanarity of multiple-planet orbits. Similarly, Sozzetti 
et al. (2001) have recently assessed the capabilities of ESA's Cornerstone Mission GAIA for the 
detection and measurement of planetary systems utilizing the v And system as a template. Our 
study is different in that, addressing the entire set of presently known multiple-planet systems, it ex- 
tends the analysis beyond the favorable cases (well-spaced, well-sampled orbits, high "astrometric" 
signal-to-noise ratios) studied by Sozzetti et al. (2001), and is specifically tailored to SIM. 

This second paper is organized as follows. In the second Section we briefly describe the setup for 
the simulation of SIM sample narrow-angle campaigns, the statistical tools implemented for planet 
detection, and the algorithms for multiple orbital fits. In the third Section we present results on 
multiple-planet systems detection, multiple-planet orbit reconstruction, and coplanarity analyses, 
and compare our results to those of previous studies. Finally, in the fourth Section we summarize 
our results and discuss our findings in the context of the present status of planet searches. 

2. Simulation Setup, Detection and Orbital Fitting Methods 

The code for the reproduction of sample observing campaigns with SIM operated in narrow- 
angle astrometric mode and the subsequent analysis of the simulated dataset has been thoroughly 
described in S02, where we defined detectability horizons and limits on distance for accurate orbital 
parameters and mass determination in the case of single planets orbiting single, nearby solar-type 
stars. In this Section we summarize its main features and working assumptions, and describe the 
modifications made in order to address the issues of multiple-planet system observations, detection, 
and measurement. 

2.1. SIM Narrow-Angle Astrometric Observing Scenario 

The SIM instrument (for a thorough presentation of the mission concept see for example 
Danner & Unwin 1999) executes pointed observations at arbitrary times and orientations of the 
interferometric baseline. When operating in the regime of narrow-angle astrometry, its fundamental 
one-dimensional measurement is the relative optical path-length delay: 



A4,n = B ■ (S* - S n ) + a d 



(1) 
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which corresponds to the instantaneous angular distance between the target and its n-th reference 
star, projected onto the baseline of the science interferometer, while the two guide interferometers 
accurately monitor changes in the satellite's attitude by acquiring fringes of two bright guide stars. 
In the above formula B = Bu^ is the interferometer baseline vector of length B, S* and S n are 
the unit vectors to the target star and its n-th reference object, while ad is the single-measurement 
accuracy on each relative delay measurement. 

Throughout all our simulations we have adopted the template observing scenario outlined in 
S02: within a domain of 1° in diameter around each target we place N r = 3 reference stars, and 
execute N Q = 24 full two-dimensional observations (corresponding to a total of N m = 144 relative 
delay measurements) randomly, uniformly distributed over the nominal L = 5 years mission lifetime, 
each composed of two one-dimensional standard visits made with orthogonal orientations of the 
interferometer baseline vector. Since in our study we have utilized the present set of extra-solar 
planetary systems, orbiting very bright (V < 10 mag), nearby, solar-type stars (spectral types F-G- 
K, except for Gliese 876), we have utilized a structure of the standard visit which applies to bright 
target and reference stars (V < 11 mag), and consequently set a single-measurement error o^, = 2 
//as throughout all our analysis (see S02 for details on the adopted error model). This value for 
is in line with the presently envisaged SIM performance in its new shared-baseline configuration. 5 

The basic assumptions made in our previous work with regard to the SIM instrument and 
the systems to be investigated are unchanged in this study: 1) we assume perfect knowledge 
of the error model and attitude of the spacecraft; 2) the objects composing the local frame of 
reference, with respect to which observations of a given scientific target are made, are assumed 
astrometrically clean, i.e. sources of astrometric noise intrinsic to the source, such as flares, rotating 
spots, or marginally unstable circumstellar disks, have not been taken into consideration, nor 
have we discussed the possibility of binaries among either targets or reference stars; 3) effects 
on the targets' displacement due to perspective acceleration, changing parallax, or higher-order 
contributions from relativistic aberration and light deflection from the major solar system bodies 
have not been included. 

We have modified the description of the target motion on the celestial sphere to allow for the 
presence of multiple planets. The stellar motion is described in terms of the 5 basic astrometric 
parameters and of the gravitational perturbations produced by the orbiting planets. The difference 
between the position vector to the target S* evaluated at time t and the same quantity S^ measured 
at time t' can then be expressed as a sum of small perturbative terms: 



where n p is the number of planets orbiting the target. The gravitational effects produced by 




(2) 



i=i 



5 The results presented and discussed in the following Sections can easily be rescaled to different timing and 
numbers of full-observations and reference stars, with the simple scaling laws already derived in S02 



-8- 



the planets are assumed to be linear (i.e., as the sum of independent Keplerian motions). Such 
simplification (no mutual gravitational interactions between planets) has in general little impact 
on the significance of our analysis. In fact, assuming the orbits are coplanar, over the time-scale of 
SIM observations (5 years), variations of the orbital parameters due to secular as well as resonant 
perturbative terms can be confidently considered negligible. This assumption may not be correct 
only in the case of Gliese 876, where the two gas giant planets orbiting the star are locked in a 
2:1 resonance, and strong mutual perturbations are already evident in the data over the timescale 
of the radial velocity observations (~ 13 years). We plan to study in detail in a future work the 
effects induced by mutual gravitational perturbations between planets on the quality of multiple- 
planet orbital fits that do and do not contain analytic descriptions of the resonant and/or secular 
perturbative nature of the systems. 



2.2. Statistical Tools for Planet Detection 

We adopt two different procedures for investigating the efficiency of SIM in detecting multiple- 
planet systems. 

First of all, following the approach in S02, we employ a standard \ 2 test with confidence level 
set to 95% applied to the observation residuals. The test is performed initially against the null 
hypothesis that there is no planet, then after removal of the first planetary signature it is applied 
in succession until no further significant deviations (at the 95% confidence level) from the fitted 
model containing n p planets appear in the residuals. 

Secondly, we have upgraded our code to include a basic periodogram analysis in order to 
search for hidden periodicities in the astrometric measurements. To this aim, we have utilized the 
standard algorithm for the evaluation of the Lomb-Scargle normalized periodogram for the spectral 
analysis of unevenly sampled data (Press et al. 1992), specifically its 'fast' version by Press k, 
Rybycki (1989). The traditional Lomb-Scargle formula is equivalent to performing an unweighted 
linear least-squares fit of a sinusoid to the data, after subtraction of the mean. For a time series 
of relative delays Ad(ti), where i takes on integer values up to the total number of measurements 
N m , the normalized periodogram as a function of the test angular frequency •& is defined as: 

2 



1 



EfiKAdfa) - Ad) cos d{U ~ t c ) E-Ll(Ac^) - Ad) sm0(U - t c ) 



^A ^ rr 



Here, Ad is the mean of the data, while t c is given by 



t c =i- arctan I ^i=i — ] ( 4 ) 



1 / sin 2^ 

— arctan — ^ 

2^ I cos 20^ 



The normalization factor is by default taken as the variance of the sample. Because of the ex- 
ploratory nature of this investigation, we assign equal weight to all observations and set ad = 2 
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fias, the value of the single-measurement error on relative delays as defined in the previous Section. 
At the moment we do not discuss the merit of other possible choices for the normalization of the 
periodogram (Gilliland & Baliunas 1987; Schwarzenberg-Czerny 1996). 

Observation residuals are inspected successively after the single-star fit, and after removal of 
each planetary signature, to establish the existence or absence of further periodicities. Also in 
this case, the significance of the detection or non detection of a periodic signal in the residuals, 
i.e. the significance of a particular peak in the power spectrum must be assessed. For the 

above choice of the normalization factor, the probability that a periodogram power z is above some 
value zq is then P(z > zq) = e~ z ° (Scargle 1982; Home &i Baliunas 1986), and, if M independent 
frequencies are scanned, the probability that none gives values larger than zq is (1 — e~ z ") M . Then, 
the false-alarm probability of the null hypothesis that the data are pure noise is computed for each 
value of z as: 

P(> z ) = 1 - (1 - e~ z °) M (5) 

This expression provides an estimate of the significance of each given peak in the spectrum 
that can be identified. For sufficiently high values of zq the false-alarm probability is very small, 
and this corresponds to the detection of a highly significant periodic signal. 

The choice of M is critical in assessing the true significance of a given peak in the spectrum. In 
general, the larger the number of frequencies, the less significant is any 'bump' in the periodogram 
power. The actual value of M depends in general on the actual frequency range scanned, and on 
the number and dumpiness of the data points available for the investigation. For this reasons, 
it is difficult to derive a closed analytical form for the number of independent frequencies. For 
period searches between the average Nyquist period (?N yq ~ 2L/N m ) and the duration of the 
data set a reasonable choice is M ~ N m (Home & Baliunas 1986). However, when sampling 
very short periods, in general frequencies much higher than the average Nyquist frequency are 
reached, then the number of independent frequencies will be higher. For any particular case, an 
actual value of M can be estimated by simple Monte Carlo analysis: holding fixed the number and 
epochs of the observations, generate synthetic datasets containing only Gaussian (normal) noise 
representative of the actual errors, then find the largest value of P(> zq) and fit the resulting 
distribution for M in Eq. 5. By doing so, Marcy & Benitz (1989) and Cumming et al. (1999) derive 
M = 1.15N m and M = 1.1757V m , respectively, and this may be compared with the analytic relation 
M = -6.4 + 1.19iV m + 0.00098 by Home & Baliunas (1986). Given the illustrative purpose 
of our analysis, we rely on results from the literature and choose a typical value M = l.20N m 
throughout our simulations. 

2.3. Procedure for Multiple-Planet Orbital Fits 

The iterative method for the solution of the non-linear systems of equations of condition is 
based on a slightly modified version of the Levenberg-Marquardt algorithm, as described in S02, 
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which ensures stability of the solution. It has been upgraded to allow for the presence of multiple, 
independent planetary perturbations. 

The photocenter motion of the target is described by an analytical model in which the re- 
computed relative delay between the target and its j-th reference star (j = 1, . . . ,N r ) is in the 
form: 

Ad+j = B • (S*(A*, /?*, fj, x ,*, At/3,*, X 2 j, X 3jh Xy, e h 7], 77)) 

l=i 

-Sj(\j, fy, Hx,j, Vf3,j, TTj)) (6) 



The five astrometric parameters (two positions, two proper motion components, and parallax) 
of the target (A*, fix,*, A*/?,*, ?r*) and each of the iV r local reference objects (Xj, (3j, fixj, fJ-pj, TTj) 
are defined in ecliptic coordinates (see S02 for details), while for the seven orbital elements needed 
to describe the orbit of each of the n p planets around the target we have utilized the Thiele-Innes 
representation, in which the semi-major axis a;, the inclination i/, the longitude of pericenter uj[ 
and the position angle of the line of nodes Qi of the l-th planet are combined to form the quantities: 

-^1,1 = ai(cosu>i cos Ui — sinuji sinO; cosi/) 

-^2,1 = sin ^ cos — cos sin $7; cos i;) 

-^3,1 = ai(cos u>i sin + sinuji cos £li cos ii) 

-^4,1 = ai(— sin u>i sin Qi + cos uji cos £li cos ii) 



In this model, the observation equations are strictly nonlinear only in the orbital period 7), the 
eccentricity e;, and the epoch of pericenter passage t\. At the end of each simulation, the classic 
orbital parameters are recomputed from the Thiele-Innes elements (see for example Heintz 1978). 
The mass M Pt \ of the l-th planet is then determined by combining the values of 7} (in years), 
7r* and the semi-major axis of the stellar orbit around the barycenter of the system a+ t \ (both in 
arcseconds), all estimated from the fit to the observed delays, and the stellar mass M*, assumed 
perfectly known, in the following approximation of the mass-function formula for M P; i <C M* : 

^p,i^(^^-J for 1 = 1,..., n p (7) 

The multiple-orbit fitting procedure also determines the formal errors on each estimated pa- 
rameter x m as: 

<r(Xm) = V Cmm, (8) 

where the C mrn 's are the diagonal elements of the covariance matrix of the fit (Press et al. 1992). 
In the case of the orbital parameters recomputed from the Thiele-Innes elements and the mass of 
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each planet, formal errors are propagated utilizing the classic error propagation formula, including 
correlations between the measured parameters: 



where the sums are extended to the m f fitted parameters from which the fc-th recomputed quantity 
Zk depends. 

As already highlighted in S02, convergence of the non-linear fitting procedure and quality of 
the orbital solutions can both be significantly affected by the choice of the starting guesses. In the 
studies presented here this issue has been left aside because, focusing on the set of presently known 
planetary systems, we have utilized as starting guesses the values of the orbital parameters for the 
planets in each system published in the literature. Nevertheless, whenever data will be inspected 
to search for planetary signatures in the absence of a priori information on the actual presence of 
planets around a given target, important questions will have to be addressed and answered, related 
to how and to what extent effective starting values, i.e., leading to successful orbital solutions, will 
be identified from the data as a function of actual instrument performances, uncertainties in the 
error model, and intrinsic properties of the planetary systems. To this end, realistic global search 
strategies will have to be tested and double-blind test campaigns conducted. Work on these issues 
is in progress and will be reported in the future. 



In addressing the problem of the detection and measurement of multiple planets with SIM, 
we have utilized the entire set of existing planetary systems known to date as representative test 
cases. The results are discussed below, as follows. First, we study SIM's efficiency in multiple- 
planet detection by means of a method that combines standard x 2 statistics, periodogram analysis 
and Fourier series expansions. Next, we estimate SIM's ability to make accurate measurements of 
multiple-planet orbits and masses. Then, we gauge the boundaries of SIM's capability to determine 
the coplanarity of the orbits of planets in systems. Finally, we quantify how beneficial it would be to 
be able to combine astrometric and spectroscopic orbital elements in order to improve the accuracy 
in full-orbit reconstruction, planet mass determination, and coplanarity measurements. All results 
reported here have been obtained assuming end-of-mission analysis of the simulated datasets. We 
do not discuss in this paper modifications to the strategy for the analysis of SIM observations at 
earlier stages of the mission, when the instrument operations are still ongoing. 




(9) 



3. 



Results 
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3.1. Detectability of Multiple-Planet Systems 

In this Section we focus on establishing the boundaries of reliable multiple-planet system 
detection with SIM for a wide range of orbital arrangements, including low values of the astrometric 
signal-to- noise ratio a/ ad and orbital periods exceeding the mission lifetime. To this end, we have 
investigated the performances of various detection methods (x 2 test, Fourier series, periodogram 
analysis) in the different regimes. Our study is primarily illustrative, and clearly does not cover 
all possible arrangements of planetary systems, such as configurations with extreme orbital shapes 
and orientations, and small amplitudes of the astrometric signal (either because of small planet 
masses, or due to large distances). 



3.1.1. Probabilities of Detection 

As discussed in S02, the easiest way to quantify SIM's sensitivity to planetary perturbations 
is to parameterize it through the orbital period T and the astrometric signal-to-noise ratio a/ ad, 
where ad is the single-measurement error on each relative delay measurement, while the astrometric 
signature a is the apparent amplitude of the gravitational perturbation induced on the observed 
star by a single companion, defined as: 

" = 2r£. ( 10 ) 

where M p , are the masses of the planet and star respectively (in solar masses), a p the semi- 
major axis of the planetary orbit (in AU), D the distance of the system from the observer (in pc). In 
the case of multiple, independent Keplerian perturbations, the orbital motion of the star about the 
barycenter will contain components reflecting the orbital periods of each of the companions, with 
the peak amplitude corresponding to the object with the larger value of the product M p x a p (for a 
given distance of the system and mass of the central star). It must be noted that Eq. 10 constitutes 
only an upper limit to the actual magnitude of the measured perturbation, without any projection 
or eccentricity effect. If the plane of a circular orbit lies in the plane of the sky, the two projections 
of the perturbation will be equal, with amplitudes given by Eq. 10. If the orbital plane is inclined 
with respect to the line of sight, there will always be some projection of the orbit on the plane of 
the sky which has the amplitude given by Eq. 10, but the amplitude at other positions angles will 
be less. If the orbit is eccentric, there could be cases where the orientation of the orbital plane 
with respect to the observer's line of sight gives rise to significantly smaller apparent amplitudes 
than indicated in Eq. 10: in the worst-case scenario, when the line of apsides is aligned with the 
line of sight the maximum measured amplitude corresponds to the orbital semi-minor axis, which 
is still equal to > 0.95 a p for moderately eccentric orbits (e < 0.3), but can drop below 0.5 a p for 
high eccentricities (e > 0.8). For very eccentric orbits, the projected amplitude of the stellar orbital 
motion can thus be up to a factor 2 or more smaller than the value calculated with Eq. 10. The 
variety of orbital arrangements reproduced by the set of planetary systems known to date does not 
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present any such case, as the eccentricity of the orbits is usually moderate (e < 0.55). The only 
exception is HD 160691 b (e = 0.8), but its astrometric signature {a = 139/ sin i //as) is sufficiently 
high not to constitute a problem in terms of ability to detect it or to determine its orbit, even in 
the worst-case scenario. 

In order to tackle the issues discussed in the previous Sections, we ran simulations of SIM 
narrow angle astrometric campaigns dedicated to the observations of the 11 stars known to harbor 
multiple-planet systems, utilizing the observing scenario sketched in Section 2.1. In particular, we 
generated sets of 500 planetary systems on the celestial sphere. The values we utilized for the 
stellar parameters and orbital elements were taken from various online catalogs in September 2002, 
as summarized in Table l. 6 The analysis carried out and results presented in this and the following 
Sections are based on a total of 990 000 simulated planetary systems. 

Two orbital elements cannot be determined by radial velocity measurements, namely i and Q. 
We simulate systems with random uniform distributions for the lines of nodes (0° < Q < 180°) and 
express detection probabilities as a function of the orbital inclination, constrained to vary in the 
range 1° < i < 90°. Planet masses are scaled with the co-secant of the inclination angle. 

An important point is the choice of the relative inclination angle between pairs of components 
in multiple systems. The relative inclination i re \ of two orbits is defined as the angle between the 
two orbital planes, and is given by the formula (see for example Kells et al. 1942): 

cosi rc i = cos i in cos i ou t + sini in sini out cos(0 out - Oin), (11) 

where i- m and i ou t, f2i n and fi ut are the inclinations and lines of nodes of the inner and outer planet, 
respectively. We choose the orbits of the multiple-planet systems shown in Table 1 to be perfectly 
coplanar, i.e. i and f2 are exactly identical for all the components in each systems. This choice has 
been adopted everywhere, unless otherwise noted. 

There are three important questions related to the use of x 2 statistics as a tool for multiple- 
planet detection that we aimed to answer: 

a) does detection probability depend on the shape and orientation of the orbits? 7 

b) is detection of low signal-to-noise astrometric perturbations (the detection limit intrinsic to 
the x 2 test is reached for a /ad — > 1) influenced by the previous removal of larger astrometric 
signals? 



For an up-to-date list of extra-solar multiple-planet systems and their orbital elements as they are known at the 
time of writing see for example Fischer et al. (2003) 

7 In previous works (Lattanzi et al. 2000; S02) it was found that, for the case of single planets, averaging over 
orbital eccentricity, the astrometric detection probability based on the y 2 test of the null hypothesis that the star is 
single did not depend upon the apparent orientation of the orbit 
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3) what is the behavior of false detections (once all planetary perturbations have been detected 
and removed) as a function of a broad range of orbital arrangements and astrometric signals? 

Figure 1 shows, for each planet in each of the 11 systems considered, the results of the x 2 test 
expressed as a function of the orbital inclination i. For each given system, the test was applied to 
the residuals to a single-star model, to a model containing the previously detected companion(s), 
and finally after the removal of all the Keplerian motions. In each system, the dominant signature 
is always easily detectable, regardless of the inclination of the orbital plane. Even long-period plan- 
ets such as 55 Cnc d, due to its large mass, are unambiguously discovered. When the companion 
producing the larger signal is removed, the remaining components, due to the usually large astro- 
metric perturbation induced on its parent star (even for edge-on configurations), are also detected 
without ambiguities, independent of the value of i, with the exception of the inner planet around 
HD 38529, the innermost planet orbiting v And, and the middle planet in the 55 Cnc system. 8 For 
the first two objects, in fact, in the edge-on configuration, the amplitude of the astrometric pertur- 
bation is of order of cr^, thus the x 2 test fails to recognize their presence about sixty percent of the 
time. 55 Cnc c has a minimum signature a ~ 4 //as, yet its detection probability in the edge-on 
configuration is a little below the 95% detection threshold due to the fact that its minimum mass 
also corresponds to the configuration in which only one dimension of the orbit is measurable by 
astrometry. As i — > 0° , towards a face-on configuration, their projected masses increase, and in all 
cases detection probability reaches 100% for a value of sini such that the correspondent astrometric 
signature a ~ 2a<i- This result confirms what had been found in our previous works (Casertano et 
al. 1996; Casertano & Sozzetti 1999; Lattanzi et al. 2000; Sozzetti et al. 2001; S02). In particular, 
we find that, as long as the orbital period does not exceed the mission lifetime, detectability of 
astrometric perturbations close to the detection limit is not hampered by the presence of one or 
more companions producing larger signals. After the larger perturbations have been fitted and 
removed, a signal-to-noise ratio a/a^ ~ 2 is sufficient for secure detection (at the 95% confidence 
level), just as if the planet was the only orbiting companion. 

The last point addressed in Figure 1 concerns the fraction of false detections after all the 
planets have been discovered and removed. Due to the confidence level adopted for the x 2 test, 
we expect ~ 5% of false positives if the subtraction of all planetary signatures has been successful. 
On the other hand, we may expect in principle that orbits with long or very short periods and 
producing low signal-to-noise ratios might not be accurately fitted and removed, with the result of 
larger spurious residuals, that could translate in a larger amount of false detections. As shown in 
Figure 1, the rate of false positives is roughly independent on the value of i, and comparable for all 
the systems investigated. Again, three exceptions are present, 55 Cnc c, HD 38529 b, and v And 
b. In these three cases, the signal to be fitted approaches the detection limit as we move towards 



8 We recall that the astrometric signature a oc T 2 / 3 , while the radial velocity amplitude K oc T -1 / 3 , thus the Hot 
Jupiters on very short-period orbits, the first to be detected by spectroscopy, are most difficult to detect with the 
astrometric method 
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i = 90°. The quality of the fits when a ~ ad is poor, and this degradation translates into a number 
of false detections up to a factor 2 larger than the expected 5%, with a more pronounced trend 
as we move close to edge-on orbits. Instead, whenever all the planets in a given system produce 
a signal a > 2 ad (even in the minimum mass configuration) , then the number of false detections 
drops even below the expected threshold. We find that the average false detection rate, including 
all systems, is = 3.2%, well within the predicted 5% of false positives in case of satisfactory 
removal of all planetary perturbations. 

3.1.2. Periodogram Analysis 

The x 2 test only allows us to draw general conclusions on the need to reject or keep an a priori 
theoretical model describing the observations, but provides no clues to the actual nature of the 
residuals, if the assumed model is proven to be not consistent with the data. In order to provide a 
direct identification of the nature of the observation residuals and determine the true significance of 
planet detection/non-detection, a standard procedure is that of inspecting the data by means of a 
periodogram analysis, to search for evidence of the presence of a clear periodicity in the time-series. 

At the end of each \ 2 test; we perform a periodogram analysis (utilizing the tools described 
in Section 2.2) to assess the presence or absence of further periodic signals. This test is applied to 
a set of observations that contain an a priori known set of planetary signatures, but nonetheless 
it has proven a useful exercise in order to gauge the ranges of applicability of this procedure to 
future, actual observations. 

Table 2 shows a set of results pertaining to simulations of the 11 planetary systems assuming 
an average value of the common inclination angle (i = 45°). The crucial points to be discussed are 
related to how well the period of the signal can be identified, and to how strong the signal must 
be with respect to the single-measurement precision, in order to give a peak z p in the periodogram 
power allowing for a significant detection, i.e. a small value of the false-alarm probability. For 
the purpose of this analysis, we set a threshold P(z p ) < 10~ 3 , which implies, given the assumed 
number of independent frequencies discussed in Section 2.2, a value of z p > 12.0 to ensure a 
significant detection of a periodic signal. 

From inspection of Table 2, the first conclusion that can be drawn is that even in the case 
of an astrometric signal-to-noise ratio ajo^ ~ 1, as long as the orbital period is shorter than the 
mission duration, a significant peak can be found in the spectrum. This is the case for example 
for HD 38529 b, that exhibits a/ad = 1-23 (for i = 45°), and while xt approaches 1 (in fact, from 
Figure 1 its presence goes undetected about 15% of the time), the false alarm probability is well 
below the threshold, and the actual peak in the spectrum matches the true orbital period to within 
~ 1%. A graphical interpretation of the results for this system is given in Figure 2. 

Secondly, three main period ranges can be identified, over which the efficacy of the periodogram 
analysis is significantly different. First of all, the orbital period may happen to be short enough 
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that the spectrum may peak at a value of T which departs significantly from the true one. This is 
for example the case of v And b, for which (see Figure 1) the probability of detection is close to 
100% (in fact, a/ad = 1-64, at i = 45°), but the period is more than a factor 3 shorter than that of 
HD 38529 b, and the periodogram analysis fails to recognize the correct value of T, identifying a 
significant peak in the spectrum at a period value which is about an order of magnitude off the true 
one. As shown in Figure 3, bottom-left panel, only the second highest peak is close to the actual 
value of the period of v And b, but still about a factor 1.5 off. The failure of the periodogram 
analysis at very short periods can be dealt with at least in part in two ways. First of all, this problem 
may be corrected by increasing the number of observations. From the results in Table 2, in fact, 
we see how, in the case of randomly spaced data, it is possible to retrieve meaningful information 
on periods much shorter than the 'average' sampling rate in the case of evenly spaced data. For 24 
two-dimensional observations distributed over 5 years, the average sampling rate would be ~ 0.2 yr, 
or 72 days. Due to the fact that some points are spaced much closer than that, the periodogram can 
effectively identify signals with periods as short as 15 days. Increasing the number of observations 
would obviously increase the ability to recover even shorter periods. Improvements in this direction 
can be obtained also by extending the range of test frequencies to values that are much greater 
than the average Nyquist frequency, and by over-sampling the spanned frequency range, in order to 
obtain information on frequency bins that are much smaller than the typically lowest independent 
frequency, that corresponds to the inverse of the timespan of the input data (i.e., the frequency 
such that the data can include one complete cycle). In our present analysis we have not quantified 
the effect of these different choices. 

Furthermore, the results summarized in Table 2 tell us that well-sampled periods up to roughly 
the mission lifetime constitute 'easy' cases, for which the correct period can be identified with 
high confidence and accuracy (typically better than 7-8%), and this independently on the actual 
magnitude of the astrometric signature, as far as it lies above the single-measurement error. Figure 4 
provides a graphical example of the results in Table 2 for the 'easy' case of HD 82943. 

Finally, when T > 5 yr we enter a regime in which the classic periodogram is intrinsically 
doomed to fail. Increasing the number of observations (without increasing the mission length) 
or manipulating the frequency range and the way it is sampled does not improve the results. In 
fact, several authors (e.g., Black & Scargle 1982; Walker et al. 1995; Cumming et al. 1999) in the 
past highlighted that the periodogram power at periods exceeding the timespan of the data can 
be significantly reduced in the traditional periodogram formula. From the results in Table 2 we 
note how, in all the cases where the actual orbital period of the planet is greater than 5 years, 
the periodogram still identifies a significant periodic signal, but the peak systematically appears at 
a period of order of the total timespan of the data. Figure 5 provides a graphical representation 
of the data in Table 2 for the classic case of 55 Cnc, which harbors the planet with the longest 
period discovered so far. The situation is somewhat different from the spectroscopic case,where as 
T lengthens the radial velocity amplitude decreases, and indeed the peaks in the power spectrum 
relative to long period signals are significantly reduced. For astrometry, the amplitude of the 
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signal increases with increasing T, and in the sample of the presently known extra-solar planetary 
systems all companions are massive (as they are easier to detect by spectroscopy). The astrometric 
signal-to- noise ratio a jod 3> 1, and so, even if only a fraction of the orbit is sampled during the 
5 years of simulated mission lifetime, the acceleration in the residuals is easily disentangled from 
the proper motion and a significant peak in the spectrum can also be identified. Nevertheless, the 
exact period cannot be recovered in a periodogram analysis. This systematic trend of a significant 
peak detected at T ~ L whenever the true period is greater than the timespan of the data is 
likely to become less evident for sufficiently small-mass planets (or, equivalently, sufficiently distant 
systems) on sufficiently long-period orbits. In fact, a reduction in the measured amplitude of long- 
period signals may occur either in the case of a/ ad — > 1 or when T is so long that secular effects 
such as proper motion begin overlapping efficiently to the astrometric signature. We have not 
investigated these issues here, as the main focus of this paper is on the actual characteristics of the 
multiple-planet systems discovered so far, and how they can be detected and measured with SIM. 



3.1.3. Search for Periodicities when T > L 

As we have seen, there exist regimes of orbital periods (shorter than the SIM mission life- 
time) for which the Lomb-Scargle periodogram approach is robust, but others (comparable to or 
longer than the duration of the observations) for which the periodogram analysis fails to provide 
sensible results. One possible solution to overcome the 'impasse' was proposed in the recent past 
(Walker et al. 1995; Nelson & Angel 1998; Cumming et al. 1999). It consists of utilizing the zero 
point of the sinusoid in the fit as a free parameter, i.e. letting the mean of the data to float during 
the fit. This is contrary to what happens in the classic periodogram formula, in which the mean 
(taken as a proxy for the zero point) is subtracted from the data, and thus assumed to be per- 
fectly known. When applied to radial velocity data, the floating mean periodogram approach has 
been shown to be successful at confirming existing companions and characterizing variations in the 
observation residuals due to long-period companions with orbits still awaiting phase closure (e.g., 
Fischer et al. 2001). We choose a different approach, with the aim of illustrating the potential of 
a strategy for period searches that tiles the classic periodogram approach to the Fourier analysis 
of periodic signals. When combined with the standard % 2 statistics, this multi-method approach 
to planet detection provides an operational framework that may also be applied more in general 
to datasets from other astrometric observatories such as the upcoming ESA Cornerstone Mission 
GAIA (Perryman et al. 2001). 

From the theory of Fourier analysis we know that a periodic, but non-sinusoidal, motion can be 
represented by a series containing a sinusoidal term at the fundamental frequency (corresponding 
to the actual period of the motion) plus higher order harmonics at exact multiple integers of this 
frequency, which carry on the information about the magnitude of the distortion from a purely 
sinusoidal motion due to the eccentricity (see for example Monet 1983). In the tangent plane to 
the celestial sphere, for the j-th planet in a system, the two components of its projected orbital 
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motion can then be written as: 

oo oo 

= ^2 a kJ cos k-fjt + ^b k j sinkjjt 

k=0 k=0 

oo oo 

Uj, tan = c k j cos k-fjt + ^ d k j sin k-fjt 

k=0 k=0 

where j = 1, . . . , n p , and 7j = 2n/Tj. The series is quickly convergent except for large eccentricities. 
For moderately eccentric orbits, say e < 0.5, three terms (fundamental, first and second harmonic) 
are sufficient to cover more than 95% of the signal (Jensen & Ulrych 1973; Konacki et al. 2002). 
For the purpose of our analysis, we have tested the ability of the above model to detect a significant 
long-period signal in the case of the five multiple systems (HD 37124, HD 38529, HD 74156, 55 
Cnc, and 47 UMa) containing one planet on a moderately eccentric (e < 0.40), long-period orbit 
(T > 5 yr). We have limited the Fourier expansion to the second harmonic term (k max = 3), and 
performed a linear fit (the unknowns being the parameters a k , b k , c k , d k , for k = 0, ... , k max ) to 
the simulated SIM astrometric data relative to the systems described above, as a function of a 
large number of trial frequencies (4000) corresponding to the period range 0.01-20 yr. Under the 
assumption of coplanarity of the orbits, then in the case of 47 UMa the inner planet has a larger 
signature, thus the outer planet signal has been searched for utilizing a model in which the Fourier 
expansion is superposed to a fully Keplerian orbit for the inner planet. The results are summarized 
in Figure 6. 

The two upper, the two middle, and the lower left panels show the behavior of the reduced 
chi-square xt as a function of the dense grid of trial periods for the abovementioned systems, again 
setting an average value of the inclination i = 45°. As opposed to the periodogram analysis, the 
Fourier series approach is capable of identifying the correct period of the outer planet, corresponding 
to the minimum of Xui to better than 4% in all the five cases. The periodicity with the larger 
amplitude can be correctly retrieved with high accuracy independently of the underlying presence of 
the second, inner companion, which is clearly evident from the fact that xt still departs significantly 
from unity even at its minimum. There are two exceptions, i.e. 47 UMa c, beyond which at present 
no other planetary companions have been discovered, and HD 38529 b, for which xt appears to 
be consistent with the outer planet being the only one in the system due to the very low signature 
induced by the inner planet. Nonetheless, as we have seen in the previous Section, the signal 
from the inner companion can still be identified in the periodogram power spectrum with high 
confidence, even in a case in which the x 2 statistics approach is more sensitive to failure. These 
results underline the importance of utilizing a range of different tools for planet detection, which 
can be combined in order to improve significantly the efficiency with which planetary signatures 
can be revealed even at low values of the astrometric signal-to-noise ratio a/ad- 

Finally, we note how in Figure 6 the minima of xl for HD 37124 c, 55 Cnc d, and 47 UMa c are 
broader than those for HD 74156 c and HD 38529 d. This is primarily due to the combined effect of 
the magnitude of the astrometric signature induced on the parent star and the fraction of the orbit 
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covered by the observations. For comparison, the lower right panel shows the results of a simulation 
in which the same fitting procedure was applied to our own solar system (nine independent, coplanar 
Keplerian orbits), placed at different distances from the observer. The dominant signature is the 
one associated with Jupiter, and its period is easily identified with an accuracy of 5-6%. As we 
move farther away from the observer, the value of the minimum of xti which is initially located 
at ^Q(min) = 1.26, indicating the presence of a second underlying signal (from Saturn, the second 
heaviest planet in our system and the one producing the second largest signature), decreases, and 
eventually at 10 pc it is perfectly consistent with Jupiter being the only planet around the Sun. 
At the same time, while the presence of the second largest planet becomes undetectable, also the 
signature from Jupiter decreases and correspondingly, for the same fraction of the orbit covered by 
the observations, the slope of the xl curve at long periods becomes increasingly less steep. 

In the perspective of the optimization of a global search strategy for the determination of 
the best configuration of starting values for the orbital parameters in the fully Keplerian fit for 
newly discovered planetary mass objects, a precise guess for the actual value of the orbital period is 
needed, in order to try to minimize the dimensions of the parametric space to be searched. Suppose 
in fact we select as plausible preliminary values of T those whose xt difference with the minimum 
Xumin ^ s l ess t nan 15% of the latter value, as done by Guirado et al. (1997) in their astrometric 
orbit determination of a low-mass companion around the star AB Doradus. This would mean that, 
for the cases shown in Figure 6, assuming no spectroscopic orbits are available, plausible starting 
guesses for the orbital period in a subsequent global least-squares fit would lie in the ranges reported 
in Table 3. As it is clear, whenever we have a configuration approaching a regime in which xt nas 
a very broad minimum (HD 37124 c, 47 UMa c, 55 Cnc d, Jupiter at 10 pc), this corresponds to a 
range of plausible starting guesses for the period covering several years, and it is likely this would 
cause accurate orbit determination and mass measurements to become a more difficult task. 

Clearly, the results shown in these Sections constitute only a very preliminary investigation 
of the potential of the approach to the problem of real-life detection of newly discovered plane- 
tary systems with a method that combines x 2 statistics, periodogram and Fourier analysis of SIM 
astrometric observation residuals. Due to the illustrative purpose of our work, some important 
topics have been left aside, such as estimates of the degradation in SIM's ability to detect plan- 
ets with low values of a/ ad on long-period orbits due to significant correlation with the target's 
proper motion. Also, some relevant technical issues remain open, such as the possibility to not only 
provide robust estimates of the reliability of the adopted model, through additional statistical hy- 
pothesis testing (e.g., Bevington & Robinson 2003), but also assess the reliability of the estimates 
themselves, through further testing and placing statistical constraints on the post-fit observation 
residuals and their covariance matrices, in order to select the optimal least-squares solution (e.g., 
Bard 1974; Bernstein 1997). Nevertheless, within the limited scope of our study, we have been able 
to provide quantitative examples which are useful for understanding that it will be very beneficial 
if a variety of numerical and statistical tools will be implemented to maximize the robustness of 
the method aimed at verifying the true significance and reliability of a detection. 



-20- 



3.2. Multiple-Planet Orbits Reconstruction 

Once a reliable detection of a system of planets has been established, a number of other 
important issues must be addressed, all primarily related to how well it is possible to measure the 
orbital parameters and masses of each planet. 

As soon as the number of planets in a system is greater than one, the parameter space to 
investigate becomes immediately very large. In this illustrative study, we have parameterized 
the ability of SIM to measure multiple-planet systems in terms of their astrometric signatures, 
orbital periods, eccentricities, and inclinations of the orbital planes. We have then investigated 
the potential of SIM to make meaningful coplanarity measurements by gauging how accurately the 
relative inclination of pairs of planetary orbits may be determined. Finally, we have quantified the 
improvement in the accuracy on the determination of orbital parameters and masses that maybe 
reached when both astrometric and spectroscopic data are available. 

3.2.1. Accuracy in Orbit and Mass Determination 

The first relevant question to be addressed is: what is the impact on the accuracy in orbital 
parameters and mass determination when the number of planets in a system is larger than one, 
and the number of parameters to be fitted correspondingly increases? 

Table 4 shows, for each planet in the 11 systems under consideration, the ratios of rms errors 
a m for the most relevant parameters (M p , a, T, e, and i) derived from a multiple-planet fit to the 
same quantities a s obtained from a single-planet solution (in simulations where only one planet was 
generated around the observed parent star). The values are averaged over the inclination angle. 

In general, for planets with periods T < L, the two-planet fit (or three-planet fit, as in the 
case of v And and 55 Cnc) degrades the orbital elements and mass estimates by ~ 30 — 40 % 
with respect to the single-planet fit, with typical fluctuations of order of 5 — 10 % among different 
parameters. From closer inspection of Table 4, some particularly interesting cases stand out. 

First of all, there is a clear trend for larger values of the ratio cr m /a s for a given parameter, 
as the magnitude of the astrometric signature decreases. For example, HD 82943 b and HD 82943 
c, with signatures of only 65/ sini //as and 22/ s'mi //as respectively (see Table 1), have values of 
the rms errors typically about a factor 2 or more larger than the respective single-planet 
opposed for example to the case of the Gliese 876 system, in which the two planets, producing 
larger signals, are measured together only about 20% worse than if they were the only orbiting 
bodies. 

The second relevant feature arising from the results of Table 4 concerns those systems in which 
the orbit of the outer planet is not fully sampled during the 5 years of simulated SIM observation. 
In three such cases (HD 38529 c, HD 37124 c, and HD 168443 c), the ratios o m ja s < 1 for almost 
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all parameters shown in Table 4. The reason for this is to be found in the intrinsic superiority of 
the model containing more than one planet, when it is applied to configurations where the outer 
planet has T > L and produces a signature much larger than the inner one. In fact, the larger 
number of free parameters (with respect to the single-planet case) is compensated by the superior 
model, and the net result is that in the least squares solution a smaller fraction of the outer planet's 
signal is identified as proper motion, with the consequence of a more precise determination of the 
outer planets' orbital elements and masses as opposed to the single-planet solution (up to about 
a factor 2-3 better). This behavior is not very prominent for 55 Cnc d, due to its much longer 
period (T ~ 3L), which translates in such poor sampling of the orbit that the larger signature is 
not enough for the orbital solution to improve significantly. 

Finally, the case of the 47 UMa system stands out. Assuming coplanarity of the two planets' 
orbits, then the astrometric signature of the outer one is smaller. The inner planet's orbit is sampled 
completely only once during the timespan of SIM observations, and the outer planet's period is 1.42 
times longer. Even utilizing the spectroscopic orbital elements as starting guesses for the iterative 
least square procedure, the net result is that a significant portion of the inner planet's signal is 
assigned to the outer one, and the accuracy in the determination of semi-major axis, period, and 
mass is degraded by over an order of magnitude with respect to the single-planet solution. For 
what concerns the uncertainties on eccentricity and inclination angle, the degradation effects due to 
this particular configuration are only marginal. However, the orbital solution for the outer planet 
does not benefit from this redistribution of the inner planet's signal, and its orbital parameters and 
mass are also more poorly determined than in the single-planet case. 

In Figure 7 and 8 we provide a detailed visualization of the results on orbit reconstruction 
and planet mass determination for the three-planet system 55 Cnc and the two-planet system HD 
12661, which are representative of the variety of configurations addressed in this study, i.e. well- 
sampled, well-spaced orbits, or poorly sampled due to either very short or very long periods, high 
values of the astrometric signal-to- noise ratio as well as values of a/a^ approaching 1. In the two 
Figures, the fractional deviation of the measured from the true value for the five relevant parameters 
discussed above is plotted against the inclination angle, both in the case of single- and multiple- 
planet realizations. Four of the five parameters (M p , a, e, and i) are primarily sensitive to the value 
of a/ ad- The larger the astrometric signal-to- noise ratio, the more accurate their measured value, 
independently on the orbital period. So for example in the 55 Cnc system the outermost planet 
has these parameters best measured (even though only one third of its orbit is sampled) , while the 
middle planet's parameters are the least accurately determined. An important exception in this case 
is the orbital eccentricity of the innermost planet: its orbit is tidally circularized (e = 0.02), and 
although its astrometric signature is a factor 2 larger than the middle planet, its period is so short 
that the inadequate sampling causes the value of e to be more poorly constrained as compared 
to the one of the middle planet. On the contrary, the accuracy with which T itself is retrieved 
primarily depends on how well the orbit is sampled (center left panels in Figures 7 and 8). Thus, 
wider orbits have the period less accurately measured, while the accuracy on the determination of 
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the other parameters is improved. 

As the inclination % decreases, we move from a perfectly edge-on to almost face-on configuration, 
i.e. from a situation in which only one dimension is actually measured to one in which the SIM 
measurements provide full two-dimensional information. At the same time, the true mass inferred 
for each planet from the radial velocity data increases, therefore the astrometric signature is larger. 
These effects combined translate into increasingly more accurate measurements of the planet's 
mass, semi-major axis, eccentricity, and period, as i — > 0°. At the same time, the observations 
are less sensitive to the inclination itself, and its fractional measurement error increases (lower left 
panels in both Figure 7 and 8). 

The present set of extra-solar planetary systems covers a relatively wide range of shapes and 
sizes of the orbits, and planetary masses. We can then attempt to generalize the results presented 
here to characterize the variety of configurations of potential planetary systems in the solar neigh- 
borhood for which SIM could provide accurate measurements, for the given template observing 
strategy outlined in Section 2.1. In the three-dimensional parameter space defined by astrometric 
signal-to-noise ratio a/ad, period T, and eccentricity e we have identified the loci for the measure- 
ment of a given orbital parameter or mass of a planet in a multiple system to a specified level of 
accuracy. Figure 9 shows iso-accuracy contours for M p , a, T, e, i, and f2 (expressed in percentages 
of the true value) in the e — T plane. The results are averaged over the inclination angle. The 
general indication is that the correlation between these two parameters is relatively limited. When 
the sampling is poor due to very short-period orbits, the degradation in accuracy is common to all 
parameters, and regardless of the orbital eccentricity. For well-sampled orbits, the parameters can 
be measured better, with varying degrees of improved accuracy. Again, the results are essentially 
insensitive to e. In fact, massive planets on relatively long-period orbits such as 55 Cnc d produce 
astrometric signals large enough that the fraction of the orbit covered by the observations still 
allows for accurate measurements of the orbital elements, regardless of the value of the eccentric- 
ity. An important exception is constituted by the accuracy on the eccentricity itself (central right 
panel), which appears to degrade for periods exceeding the 5-year SIM mission duration, unless 
e > 0.1. 9 Finally, it must be noted how the accuracy in the determination of the orbital period 
(central left panel) shows a behavior opposite to that of the other parameters, in agreement with 
the fact that how well the periodicity of the astrometric signal can be measured depends upon the 
number of times an orbit is fully sampled. The same mostly uncorrelated behavior is evident in 
the plane defined by signal-to-noise ratio and eccentricity. 

Figure 10 shows the same results, but this time in the plane defined by a /ad and T. Here we 
appreciate a clear correlation between these two crucial parameters, as it is shown by the trend 



9 When the orbital period is sufficiently long or the mass of the planet sufficiently low, instead, e is likely to play 
an important role in the degradation of the reachable accuracy on other parameters as well, as eccentric orbits not 
sampled in correspondence of the pericenter passage may not be properly measured by the fitting algorithms due to 
a lack of curvature sufficient to disentangle the periodic motion from the stellar proper motion 
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for lower accuracy on M p , a, e, i, and towards low values of a/ad and short periods, and for 
improved determination of orbital elements and masses for large values of a/ ad and longer periods. 
Again, the orbital period itself behaves in an opposite fashion. 

It is possible to provide an analytic representation for the dependence of the fractional errors 
on the estimated orbital elements and masses as a function of the three quantities a/ ad, T, and 
e. In order to do so, we have performed a fit to the data (i.e., the fractional uncertainties a p on a 
given parameter, where the subscript p can be equal to M p , a, T, e, i, or Q) with a model of the 
form: 

a p (%) = m x (a/ad f 1 x T m x e" 3 (12) 

The results are summarized in Table 5, and are averaged over the inclination angle. Again, similarly 
to what we saw in Figures 9 and 10, we realize how, for example, the fractional errors on M p , a, 
and i are dominated by the behavior of a/ ad and T, while the dependence on e in such cases is 
relatively weak. For a e , instead, the dependence on the eccentricity itself is much stronger. Finally, 
in the case of <tt, the exponent of a /ad, that would have the fractional error on the orbital period 
decrease as the astrometric signal-to-noise ratio increases, is largely compensated by the opposite, 
and stronger, dependence on T itself, with the result that, as T increases, <tt grows larger as well. 
It is important to stress the limits of applicability of such empirical multi-dimensional power-laws. 
In fact, these formulas do not include a variety of effects, such as the degradation induced by 
the increasingly poor sampling at very short or very long periods. Due to the former effect, for 
example, the best-fit formula for <tt systematically underestimates the errors on the orbital period 
for T < 0.1 yr. In general, the agreement between the power-laws shown in Table 5 and the data 
is good when the values of a p are larger than 1%, while discrepancy factors of a few typically arise 
when the fractional errors on the parameters become very small. Finally, we expect the functional 
dependence of the fractional errors on a/ad, T, and e, as well as the values of the zero-points 
7/o, to be somewhat different when a range of possible observing strategies and different values of 
the measurement errors were to be considered. In our exploratory work we have not investigated 
further such issues. 

Taking into account the results reported in Tables 4 and 5 and those summarized in Figures 9 
and 10, we can attempt to draw some general conclusions as for what concerns the ability of SIM 
to measure the most relevant parameters in multiple-planet systems: on average, 10% accuracy, or 
better, will be reached in the determination of the mass of planets in two or three-planet systems 
with components on periods not shorter than T ~ 0.1 yr, and producing an astrometric signal-to- 
noise ratio a /ad > 5. Similar accuracy will be attainable for systems containing massive planets 
with periods longer than the SIM mission duration (up to T ~ 3L), which produce a /ad > 100. The 
semi-major axis, orbital inclination, and lines of nodes will behave similarly, while the eccentricity 
will require a/ ad — 10 for periods T < L, and accurate measurements will become increasingly 
more difficult to achieve for periods exceeding the mission lifetime. The orbital period will be 
recovered with typical accuracies of a fraction of a percent for T < L. Even for objects with 
T ~ 3L and producing a/ad ^ 1 ; such as 55 Cnc d, 1% accuracy should be attainable. 



3.2.2. Coplanarity Measurements 



As discussed in the Introduction, the increasing number of extra-solar planetary systems has 
motivated detailed theoretical studies on their dynamical evolution and long-term stability. It is 
worth mentioning in detail the main results derived so far for some of the systems which have been 
the target of the larger theoretical efforts. For example, in the case of the v And system, neglecting 
to first order the effects of the innermost planet on the overall stability of the system, Stepinski et 
al. (2000) come to the conclusion that dynamical stability requires the orbital inclination of the 
outer two companions to be greater than % ~ 13°, otherwise the two objects would be too massive 
and gravitational interactions would disrupt the system. Furthermore, the system cannot be stable 
in the long term if relative inclinations are greater than 55°, 35°, and 10°, for i ~ 64°, i ~ 30°, 
and i ~ 15°, respectively: then, the more massive the planets, the closer to coplanarity their orbits 
have to be, for the system not to be destabilized on a short time-scale. As for what concerns 
the strongly interacting two-planet system on short-period orbits around Gliese 876, Laughlin & 
Chambers (2001), Rivera & Lissauer (2001b), and Gozdziewski et al. (2002) have found limits for 
dynamical stability on their relative inclinations that, similarly to the v And case, are a strong 
function of the unknown inclination of the orbital planes (and thus unknown masses): the relative 
inclinations of stable systems can vary in the ranges of ±15° for sini = 0.5 up to ±90° for sini = 1. 
Finally, Laughlin et al. (2002) have studied extensively the 47 UMa system, and found that, 
somewhat surprisingly, the mutual inclination of the two long-period planets orbiting 47 UMa 
must be less than 45° for the system to be stable, but the results are essentially insensitive to the 
actual planet masses. 

As the above results clearly suggest, dynamical limitations on the relative inclinations and 
masses of planetary companions in extra-solar multiple-planet systems cannot be stated very pre- 
cisely. Such estimates suffer from large uncertainties due to the typically poorly constrained orbital 
inclinations and planet masses, and undetermined position angles of the lines of nodes. More 
accurate observational data to complement the one-dimensional radial velocity measurements are 
needed before the details of the long-term evolution of multiple-planet systems can be assessed with 
a high degree of confidence. 

Figure 11 shows the estimated accuracy with which SIM would measure the coplanarity of 
the orbits of the present set of extra-solar planetary systems. After the multiple-planet fit, the 
estimated values of i and Q for all components in each system were used to derive their actual 
relative inclination (Eq. 11). In Figure 11, for each system, the relative inclination of each pair of 
planets is plotted against the (common) inclination of the orbital plane with respect to the plane 
of the sky. The corresponding error bars have been computed by propagating (Eq. 9) the formal 
expressions from the covariance matrix of the multiple-planet fit. 

In order to establish coplanarity, we must measure i rc \ ~ 0°. The first seven panels show 
how, as long as all components in a given system have minimum masses and orbital periods such 
that the astrometric signal-to-noise ratio is favorable (a/ad — 10 or greater) for any inclination 
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of the orbital plane with respect to the plane of the sky, the quasi-coplanarity of each pair of 
planetary orbits can be assessed with high accuracy. Even in the case of two systems with one of 
the planets having period exceeding the timespan of the observations (HD 37124 and 47 UMa), 
the net result is that for such well-measured systems the relative inclination can be determined 
to be i rc \ < 3°, with relative uncertainties of a few degrees. The next seven panels of Figure 11 
highlight instead how, when dealing with astrometric signatures approaching the detection limit 
(a/ad ~^ 1)) to make accurate coplanarity measurements will be a significantly more challenging 
task. In the specific cases shown, only towards quasi-face-on configurations, in which the signal 
from the smaller and/or shorter period planet becomes large enough (due to significantly larger 
projected mass), the results resemble those obtained in the favorable cases. Otherwise, the relative 
inclinations will typically be measured to be i rc \ ~ 10° or greater, with uncertainties of additional 
tens of degrees. Finally, as the bottom right panel shows, for the two outer planets of the v And 
system quasi-coplanarity can be reliably established, in agreement with the findings of Sozzetti et 
al. (2001). 

It is worth mentioning another interesting feature arising from the results presented in Fig- 
ure 11. Almost all panels show a more or less pronounced trend for increasing uncertainties in 
the determination of i re \ as we move away from the edge-on configuration, with a somewhat sharp 
change in behavior as % — ► 0°. The explanation is to be found in the details of the dependence of 
i rc \ on i; n , i out , and the difference f2 out — fii n - In fact, as i decreases, the uncertainty in the position 
angle of the line of nodes grows, because of the increasing difficulty in its correct identification (Q 
eventually becomes undefined at i = 0°). However, the uncertainty in the value of i rc \ defined in 
Eq. 11 increases correspondingly only up to the point in which, for sufficiently low values (i < 5°) of 
the inclination angles, the second term sini; n sini out cos(f] out — ^; n ) of the right-hand side becomes 
essentially negligible with respect to the first term cos i- m cos i out , irrespective of the fact that the 
uncertainty of the line of nodes keeps on growing. Thus, ultimately for quasi-face-on configurations 
an accurate knowledge of Cli n and O ou t is not required. In their recent work on extra-solar plane- 
tary system detection and measurement with GAIA, Sozzetti et al. (2001) do not consider orbital 
configurations of the v And system with inclination angles less than 5 degrees. As a consequence, 
they come to the conclusion that when i — > 0° accurately estimating the coplanarity of the outer 
two planets of the v And system would be difficult because of the increasing uncertainties on the 
measurement of for the two planets. Here instead we have shown how the relative inclination 
becomes essentially insensitive to the retrieved values of Q- m and £l ut, as we approach a perfectly 
face-on configuration. 

Finally, simulations of non-coplanar systems with f2; n / fi out and/or i- m / i out yielded similar 
results. Our findings help to reaffirm the importance of SIM high-precision position measurements 
in verifying the stability of multiple-planet system. Accurately determining the inclination angles 
and lines of nodes of multiple planetary orbits will allow in turn to derive meaningful mass and rel- 
ative inclination angle estimates, which will be used to better constrain the results from theoretical 
studies on the long-term evolution of planetary systems. 
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3.2.3. Combining Radial Velocity and Astrometric Data 

In recent works, Eisner & Kulkarni (2001a, 2001b, 2002) have utilized a semi-analytical method 
to show how planet detection efficiency would benefit from the simultaneous availability of both 
astrometric and radial velocity measurements, especially as far as long-period, edge-on orbits are 
concerned. This combined approach has indeed proved very successful when dealing with real data. 
For example, utilizing both astrometric measurements from the Fine Guidance Sensor 3 on board 
the Hubble Space Telescope and ground-based precision radial velocity data, Benedict et al. (2001) 
have improved the accuracy on the mass estimates for the M dwarf binary Wolf 1062 by a factor 
4 with respect to the same values obtained by Franz et al. (1998) using only HST astrometry. 
Very recently, the same combined method revealed itself an essential ingredient in the spectacular 
determination of the first actual mass of an extra-solar planet, the outer component in the Gliese 
876 system (Benedict et al. 2002). 

Usually, a combined astrometric + spectroscopic orbit is obtained by modeling time-series of 
data from both techniques in a simultaneous least-squares solution, with a few additional con- 
straints, the most important of which is the identity (Pourbaix & Jorissen 2000): 

asini TK\ y/l — e 2 , , 

7T* ~ 2vr x 4.7405 ' ^ ' 

which relates quantities only derived from astrometry (inclination angle, astrometric signature, 
and stellar parallax on the left-hand side) to quantities derivable from both techniques or radial 
velocity alone (orbital period and eccentricity, and radial velocity semi-amplitude of the primary 
on right-hand side). 

In our exploratory work, we have adopted a more limited approach. For each of the 11 
planetary systems under study, we have utilized the orbital elements derived from spectroscopy 
(T, e, r, and uj) as constraints in the sense that we have kept them fixed to their input values 
in the global least-squares iterative solution, and solved only for a, i, £1. Figure 12 shows the 
estimated improvement that may be obtained in the determination of planet mass, semi-major 
axis, and inclination when SIM relative astrometry of the presently known planetary systems is 
"combined" (in the abovementioned sense) with spectroscopy, as opposed to the scenario in which 
all orbital elements are solved for in the least-squares minimization procedure. The three contour 
plots identify, in the plane defined by a /ad and T, regions of increasingly larger values of the 
ratio ogI°~C between the estimated fractional uncertainties on M p , a, and i in the case of 'global' 
orbital fits in which all parameters were adjusted and the same quantities as computed after the 
'constrained' fits. The results are averaged over the inclination angle. 

The general indication is that for well-sampled orbits, with large astrometric signatures the 
constrained fit does not improve significantly the final accuracy of the results. In such cases, 
typically 1 < (Jg/cc < 2. The ability to recover accurate values of M p , a, and i will be improved 
especially when the period is larger than the timespan of the observations (and the astrometric 
signal-to- noise ratio is very large), with typical values of o-g/°~c > 6 for a and M p , while the 
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inclination is less sensitive to the constrained solution. For low values of a/ ad and short-period 
orbits the improvement on the fractional accuracy for the fitted parameters will be marginal. 

In their work on multiple-planet system detection and measurement with GAIA, Sozzetti et 
al. (2001) suggested that, for systems with configurations very close to face-on, the accuracy of the 
inclination measurement will be substantially increased by combining radial velocity and astromet- 
ric data. On the other hand, although in their work they focused mainly on combining astrometry 
and radial velocity in the case of edge-on (single-planet) orbits, Eisner and Kulkarni (2002) argue 
that the short-period accuracy of astrometry + radial velocity should be approximately independent 
of orbital inclination. In Figure 13 and 14 we show the behavior of the ratios <jq/ ac for M p (upper 
panels), a (central panels), and i (lower panels) as a function of the common inclination of the 
orbital plane for four representative systems (55 Cnc, 47 UMa, HD 12661, and HD 38529), whose 
characteristics allow us to provide quantitative answers to the above arguments. Although the 
predictions made by Sozzetti et al. (2001) and Eisner & Kulkarni (2002) are in general confirmed, 
when it comes to multiple-planet orbital fits some details become more subtle and complex. In fact, 
not only the average value of the ratios ac/ac ( as shown in Figure 12), but also their behavior as a 
function of the orbital inclination depends on the ranges of periods and astrometric signatures we 
are dealing with. In the case of 55 Cnc d, for example, with a period T ~ 3L, the accuracy on all 
fitted parameters benefits greatly from a constrained fit, which minimizes the residual covariance 
between the orbital solution and the solution for the parallax and the proper motion of the primary. 
The improvement in orbit reconstruction reaches its maximum for an almost face-on configuration. 
On the other hand, for 47 UMa c and HD 38529 c, with T ~ 1.4L and T ~ 1.2L respectively, a G /a c 
increases as i — > 0° only in the case of i itself, while the semi-major axis is almost insensitive to 
the inclination angle. In the range of periods T < L, again different parameters behave slightly 
differently, in these cases primarily due to the magnitude of a /ad- While the accuracy on orbital 
inclination is only marginally improved towards face-on configurations, the estimated semi-major 
axis (and in turn the derived planet mass) is more accurately retrieved in a constrained fit on 
an edge-on configuration, in which case astrometry loses completely the second dimension of the 
measurements. 10 For an accurate determination of the semi-major axis, the combination between 
astrometric and spectroscopic data will be effective particularly when a/aa — > 1, as in the case of 
HD 38529 b or 55 Cnc c, while for example the results for the two planets in the HD 12661 system, 
both with large signatures, are only marginally improved by the constrained fit. 



10 A noticeable exception is constituted by 47 UMa b, that exhibits significant improvement in the accuracy with 
which its semi-major axis can be determined when the configuration is face-on. This behavior is due to the particular 
orbital arrangement of the system, with the longer period planet producing the smaller astrometric signal, under the 
assumption of coplanar orbits 
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4. Summary and Conclusions 

The properties of extra-solar giant planets detected by radial velocity surveys of nearby solar- 
type stars (Marcy et al. 2002a) seem to indicate that the Solar System configuration is just one 
of the many possible outcomes of disk evolution around young stars. Indeed, the strong coupling 
between the early evolutionary processes of a star and its disk can have a significant impact on 
planet formation time-scales and the final orbital configurations after the disk dissipates. 

Over the next decade of so, a series of new instruments will come on line, which will provide 
data of great value to shed new light in the complex scenarios of the formation and evolution 
of planetary systems. Among indirect detection techniques, ground-based precision spectroscopy, 
the most successful technique so far, will be complemented by high-precision ground-based (Keck 
Interferometer, VLTI) and space-borne (SIM, GAIA) astrometry, and transit photometry from 
ground (e.g., OGLE III, STARE, STEPSS, Vulcan Camera Project) and in space (Corot, Kepler, 
Eddington). In the field of direct detection techniques, near- and far- infra-red diffraction- limited 
ground- and space-based imaging (ALMA, SIRTF, JWST) will pave the way to ambitious projects 
of coronagraphic/interferometric imaging from space (TPF, Darwin), with the long-term goal of 
directly imaging terrestrial planets in the Habitable Zone of nearby stars. 

In this paper we have completed the analysis begun in our previous works (Casertano & 
Sozzetti 1999; S02), in order to connect and relate the basic capabilities of the Space Interferom- 
etry Mission (SIM) to the properties of extra-solar planetary systems. We have utilized detailed 
end-to-end numerical simulations of sample SIM narrow-angle astrometric observing campaigns 
(Section 2.1), an improved methodology for planet detection which combines % 2 statistics, peri- 
odogram analysis, and Fourier series expansions (Sections 2.2), an upgraded analytical model that 
allows for multiple-planet orbital fits (Sections 2.3), and the set of presently known extra-solar 
planetary systems as templates. The experiments described in this paper have allowed us to quan- 
tify the limiting capability of SIM to discover systems of planets around solar-type stars in the 
solar neighborhood (Sections 3.1.1, 3.1.2, and 3.1.3), measure their orbital properties and masses 
(Sections 3.2.1 and 3.2.3), and accurately determine the coplanarity of multiple-planet orbits (Sec- 
tion 3.2.2). Our main results can be summarized as follows. 

1. Additional planets in systems have little impact on SIM ability to detect each component in a 
system, in comparison to the single-planet configurations. The inaccurate fit and subtraction 
of orbits with astrometric signal-to- noise ratio a /ad — ► 1 can on the other hand increase the 
false detection rate by up to a factor 2. The periodogram analysis adds robustness to the 
detection method when T < L, by singling out periodicities even in the case of ot/od close to 
the x 2 detection limit. For very short-period orbits, a more dense time-series of observations 
will be the obvious choice in order to overcome poor sampling. When T > L, the least squares 
technique combined with Fourier analysis can correctly identify periods as long as three times 
the timespan of the observations. This approach is arguably preferable, as opposed to the 
periodogram method, which needs modifications in the long-period regime. 
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2. Accurate measurements of multiple-planet orbits and determination of planet masses are 
only moderately affected by the presence of more than one object in a system, with typical 
degradation of 30-40% with respect to single-planet solutions. When T < L, it is possible 
to determine masses and orbital inclinations to better than 10% for systems with planets 
having periods as short as 0.1 yr, and for systems with components producing astrometric 
signals as low as ~ 5 ad, while a/ ad ~ 100 is required in order to measure with similar 
accuracy objects with periods as long as three times the mission duration. Orbital eccentricity 
typically requires larger signals for the same accuracy level, and its correct identification can 
become a non trivial task when T > L. The accuracy on estimated orbital elements improves 
significantly as we move towards face-on orbits, except for the inclination angle. 

3. Accurate coplanarity measurements are possible for systems with all components producing 
ct/&d — 10 or greater. In the case of truly coplanar systems, the relative inclination between 
pairs of planetary orbits is measured to be z re i < 3°, with uncertainties of a few degrees, 
for periods 0.1 < T < 15 yr. In systems where at least one component has a/ ad — ► 1, 
uncertainties on i re i of order of 30° — 40°, or larger, are likely to preclude a robust assessment 
of the system coplanarity. 

4. Whenever feasible, an approach that combines astrometry and radial velocity will yield sig- 
nificantly more accurate estimates of planet masses and orbital elements. The uncertainties 
on orbital elements and masses can be reduced by up to an order of magnitude, especially 
in the case of long-period orbits in face-on configurations, and for low amplitude orbits seen 
edge-on. Well-sampled, well-measured orbits (T < L, a/ ad S> 1) are only marginally affected 
by the combination of astrometric and radial velocity measurements. 

Our results reaffirm the important role future high-precision space-borne astrometric missions 
promise to play in the realm of extra-solar planets. With its unprecedented astrometric precision, 
SIM will not only complement other on-going and planned spectroscopic, astrometric, and pho- 
tometric surveys, but its position measurements will have a unique impact in the study of some 
important aspects of multiple-planet systems. By directly measuring the line of nodes and incli- 
nation angle for each component in a system, SIM will determine whether planetary orbits are 
coplanar with uncertainties of a few degrees. For instance, this will provide theory with the obser- 
vational evidence needed to address the long-term evolution issue, and draw sensible conclusions 
on the chaotic and stable/unstable dynamical behavior of multiple-planet orbits. The same SIM 
measurements will also be instrumental to help confirm or rule out one of paradigms which form 
the basis for present-day theoretical models, i.e. that the sole conceivable environment for planet 
formation are flattened circumstellar disks around young pre-main sequence stars. If multiple- 
planet orbits are found to be coplanar, this will indicate that planetary systems indeed originate 
and evolve in a way similar to our own (Lissauer 1993; Pollack et al. 1996). If large relative orbital 
inclinations were found, this would provide evidence that other systems are truly different from 
ours, and thus their present configurations should be explained in terms of either an early, chaotic 
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phase of orbital evolution or formation by another mechanism such as disk instability (Boss 2000, 
2001; Mayer et al. 2002). 
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Fig. 1. — Probability of detection as a function of orbital inclination i, for each planet in each of the 
11 multiple-planet systems considered in our study. The three planets HD 38529 b (dashed-dotted- 
dotted line), v And b (long-dashed line), and 55 Cnc c (dashed-dotted line) are athe most difficult 
to detect. After the not accurate subtraction of their low-amplitude astrometric signatures, the 
corresponding false detection rates are up to a factor 2 larger than the expected 5% 

Fig. 2. — Periodogram analysis for the HD 38529 planetary system. In the upper left panel we 
show the periodicity search after the single-star fit: a long-period signal is clearly found with high 
significance, but due to the incomplete orbital sampling, its value is about a factor 2 shorter than 
the actual period of HD 38529 c. The upper right panel highlights how, after the removal of the 
signal with the larger amplitude, the period of HD 38529 b is still accurately identified, in spite of 
the very low astrometric signature induced by the inner planet of the system on its central star. 
The lower left panel present the results of the periodogram analysis after the dual Keplerian fit to 
the observations: no significant peaks are found 

Fig. 3. — Same as Figure 2, but for the v And planetary system. The periodogram analysis correctly 
identifies the periods of the two outer planets (two upper panels), and rules out the presence of 
additional planets after subtraction of all planetary signatures (lower right panel). Due to very 
short period of the innermost planet, instead, after removal of the signatures of v And c and v And 
d the periodogram search identifies a significant peak in the spectrum which does not correspond 
to the true period of v And b (lower left panel) 

Fig. 4. — Same as Figure 2, but for the 'easy' case of the HD 82943 planetary system 

Fig. 5. — Same as Figure 2, but for the case of the 55 Cnc planetary system. The planet inducing 
the larger signature, 55 Cnc d, is identified immediately in the power spectrum (upper left panel), 
but, similarly to the case of HD 38529 c shown in Figure 2, the incomplete orbital sampling causes 
the periodogram analysis to identify an incorrect value for the orbital period of 55 Cnc d 

Fig. 6. — Reduced chi-square xl as a function of a dense grid of trial periods for the five planetary 
systems harboring one planet with a period exceeding the timespan of SIM observations (two upper, 
two middle, and lower left panels). For comparison, the lower right panel shows the results of the 
same procedure applied to our own solar system, placed at increasing distance from the observer 
(details in the text) 

Fig. 7. — Fractional errors (%) for M p , a, T, e, and i, as a function of the inclination of the orbital 
plane, in the case of the 55 Cnc planetary system. In each panel, the results for a given parameter 
are shown for all the system components, for the cases of a single-planet solution (dashed-dotted, 
dotted, and short-dashed lines for 55 Cnc b, 55 Cnc c, and 55 Cnc d, respectively) and a complete 
three-planet model (solid, long-dashed, and dashed-dotted-dotted lines for 55 Cnc b, 55 Cnc c, and 
55 Cnc d, respectively) 
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Fig. 8. — Same as Figure 7, but for the case of the HD 12661 planetary system. The dashed-dotted 
and dashed-dotted-dotted lines are relative to the single-planet solution for HD 12661 b and HD 
12661 c, respectively. The solid and long-dashed lines are relative to the dual Keplerian fit for HD 
12661 b and HD 12661 c, respectively 

Fig. 9. — Iso-accuracy contours in the e — T plane, for the most relevant orbital parameters and 
planet mass. The results in each panel have been obtained utilizing the 11 planetary systems known 
to-date. The contour regions are color-coded by the accuracy achieved in the determination of the 
given parameter, expressed as a fraction (%) of its true value 

Fig. 10. — Same as Figure 9, but with the contour regions identified in the a/ ad — T plane 

Fig. 11. — Relative inclination i TC \ between pairs of planetary orbits as a function of the common 
inclination angle with respect to the line of sight, for the entire set of presently known multiple- 
planet systems. In each panel, the corresponding uncertainties are computed utilizing the formal 
expressions from the covariance matrix of the multiple Keplerian fit 

Fig. 12. — Contours identifying regions, in the a/a^ — T plane, with equal values of the ratios 
ac/ac of the estimated fractional uncertainties on M p , a, and i in the case of 'global' fits in which 
all parameters were adjusted to the same quantities as computed after the 'constrained' fits (as 
discussed in the text), for the set of presently known planetary sytems. The results are averaged 
over the inclination angle (coplanar orbits are assumed). The contour regions are color-coded by 
the value of the ratio gq / &c f° r each parameter 

Fig. 13. — The behavior of the ratios <jg/&c for M p (upper panels), a (middle panels), and i 
(lower panels), expressed as a function of the common inclination angle, for the two representative 
multiple-planet systems 55 Cnc and 47 UMa 



Fig. 14 — 



Same as Figure 13, but for the HD 12661 and HD 38529 planetary systems 
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Table 1. Fundamental stellar characteristics, astrometric parameters and orbital elements of the 
presently known set of extra-solar planetary systems utilized in the simulations. Planet masses 
and astrometric signatures are computed as lower limits corresponding to edge-on configurations 
(sini = 1). The data have been taken utilizing various web resources containing up-to-date values 

of all parameters, as explained in the notes 
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Table 1 — Continued 
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Table 1 — Continued 



Stellar Parameters 8, 


Central Star 


Orbital elements 13 


Planet b 


Planet c Planet d 






System: HD 12661 






A (deg) 


37.70 


a (//as) 


47.43 


100.4 


(deg) 


12.26 


a (AU) 


0.82 


2.56 


fix (mas/yr) 


-161.28 


T(d) 


263.3 


1444.5 


fip (mas/yr) 


-127.78 


e 


0.35 


0.20 


7T (mas) 


26.91 


r(JD) 


2459943.7 


2459673.9 


V Magnitude 


7.43 


w (deg) 


292.6 


147.0 


Spectral Type 


G6V 


Mpsini (Mj) 


2.30 


1.56 


d (pc) 


37.16 


Q (deg) 


[0,7T] 


[0, 7T] 


M* (M ) 


1.07 


i (deg) 


[0,tt/2] 


[0,vr/2] 






System: HD 160691° 






A (deg) 


267.20 


a (//as) 


155.8 


139.2 


P (deg) 


-28.87 


a (AU) 


1.48 


2.30 


// A (mas/yr) 


-20.96 


T(d) 


637.3 


1300.0 


A*>9 (mas/yr) 


-190.61 


e 


0.31 


0.8 


7r (mas) 


65.36 


r (JD) 


2450959.0 


2451613.0 


U Magnitude 


5.15 


w (deg) 


320.0 


99.0 


Spectral Type 


G3IV/V 


Mpsin/ (Mj) 


1.74 


1.0 


d (pc) 


15.30 


Q (deg) 


[0,7T] 


[0.7T] 


M* (M ) 


1.08 


i (deg) 


[O.tt/2] 


[O.tt/2] 






System: 47 UMa 






A (deg) 


149.30 


a (//as) 


366.1 


195.5 


/? (deg) 


31.28 


a (AU) 


2.09 


3.73 


MA (mas/yr) 


-259.16 


T(d) 


1089.0 


2594.0 


(mas/yr) 


188.90 


e 


0.06 


0.10 


7r (mas) 


71.02 


r (JD) 


2453622.0 


2451363.5 



-41 - 



Table 1 — Continued 



Stellar Parameters 21 


Central Star 


b 

Orbital elements 


Planet b 


Planet c Planet d 


V Magnitude 


5.03 


to (deg) 


172.0 


127.0 


Spectral Type 


GOV 


Mpsini (Mj) 


2.54 


0.76 


d (pc) 


14.08 


Q (deg) 


[0.7T] 


[0,7T] 


M± (Ma) 


1.03 


i (deg) 


[0,7r/2l 

L 5 / J 


[0,7r/2l 

L 5 / J 






System: HD 82943 d 






A (deg) 

V 0/ 


150.22 


a (/uas) 


65.58 


22.28 


3 (deg) 


-24.79 


a (AU) 


1.16 


0.73 


(mas/vr) 


-58.14 


T (d) 


444.6 


221.6 


Lin (mas/vr) 


-164.07 


e 


0.41 


0.54 


7r (mas) 


36.42 


r (JD) 


2451620.3 


2451630.9 


y Magnitude 


6.54 


u; (cleg) 


96.0 


138.0 


Spectral Type 


GO 


Mpsini (Mj) 


1.63 


0.88 


d (pc) 


27.46 


Q (deg) 


[0,7T] 


[0,7T] 


M* (M ) 


1.05 


i (deg) 


[0,tt/2] 


[0,7r/2] 






System: HD 74156 d 






A (deg) 


131.71 


a (^as) 


6.35 


383.9 


/3 (deg) 


-12.86 


a (AU) 


0.276 


3.47 


Ma (mas/yr) 


-28.19 


T(d) 


51.61 


2300.0 


(mas/yr) 


-200.05 


e 


0.65 


0.40 


7r (mas) 


15.49 


r (JD) 


2451981.4 


2450849.0 


y Magnitude 


7.62 


w (deg) 


183.7 


240.0 


Spectral Type 


GO 


Mpsini (Mj) 


1.56 


7.5 


d (pc) 


64.56 


n (deg) 


[0,7r] 


[0,7T] 


(M ) 


1.05 


i (deg) 


[0,tt/2] 


[0,tt/2] 



-42 - 



Table 1 — Continued 



Stellar Parameters 3, 
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a Data from http://simbad.u-strasbg.fr/. Positions and proper motion components, origi- 
nally given in equatorial coordinates, have been projected to the ecliptic reference frame (see 
for example Green 1985) 

b Data from http://www.obspm.fr/encycl/encycl.html and http://exoplanets.org/ as of 
September 2002 

c The existence of the outer planet around HD 160691 is still doubtful. Its presence has been 
inferred by Jones et al. (2002), due to a clear trend in the radial velocity residuals, but the 
orbit still awaits phase closure 

d The discovery of the planetary systems around HD 74156 and HD 82943 has not yet been 
reported in refereed Journals 



-43- 



Table 2. Results from the periodogram analysis of the 11 planetary systems known to date, in 
the case of an average value of the orbital inclination angle (i = 45°). The columns display, 

respectively, the name of the system, the type of fit performed (single-star, single- or 
multiple-planet), the value of the reduced chi-square xt °f the fit, the highest peak z p in the 
periodogram power spectrum, the value of the orbital period T(z p ) corresponding to the peak in 
the spectrum, and the false-alarm probability P(z p ) of the peak 
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42.8020 


0.157193 


0.0 




2 Planets 


0.943551 


6.85225 


0.109261 


0.146701 


HD 12661 


Single Star 


67.2261 


28.1816 


3.30105 


1.23 x 10~ 7 




1 Planet 


481.436 


44.0427 


0.741820 


0.0 




2 Planets 


1.13029 


6.83758 


0.0402566 


0.148699 


HD 160691 


Single Star 


225.807 


50.5789 


1.65052 


0.0 




1 Planet 


1650.87 


56.9019 


3.27908 


0.0 




2 Planets 


0.857078 


7.96351 


0.397465 


0.0508595 


47 UMa 


Single Star 


688.689 


42.7118 


2.81892 


0.0 




1 Planet 


402.954 


51.3840 


4.95227 


1.35 x 10~ 9 




2 Planets 


1.04254 


4.81968 


0.1009260 


0.703383 
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Table 2 — Continued 



System 


Type of Fit 




Zp 


T(z p ) 




HD 82943 


Single Star 

1 Planet 

2 Planets 


29.2090 
36.8078 
1.01055 


53.8883 
43.9404 
5.77669 


1.23789 
0.618946 
0.0398152 


0.0 
0.0 
0.372220 


HD 74156 


Single Star 

1 Planet 

2 Planets 


224.126 
4.25179 
0.994640 


48.2071 
36.1520 
5.25626 


3.30105 
0.140470 
0.0330859 


0.0 
0.0 
0.543543 


HD 37124 


Single Star 

1 Planet 

2 Planets 


82.4841 
36.8113 
0.931258 


57.8993 
39.9654 
5.95232 


4.36463 
0.415679 
0.0394809 


0.0 
0.0 
0.321418 
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Table 3. Ranges of plausible starting guesses for the orbital period of the 5 planets with T > 5 
yr in the presently known extrasolar multiple-planet systems. For reference, the same results are 
reported for Jupiter in simulations of the solar system as 'seen' by SIM at increasing distances (in 
pc) from the observer. The results are derived in the context of the Fourier analysis approach 
discussed in the text. Preliminary values of T are considered acceptable if the corresponding xt 
does not differ from the minimum xtmin by more than ±15% of the latter value 



Planet Plausible xt Ranges Plausible T Ranges 



HD 37124 c 


28.62 


< 


Xp 


. < 


38.72 


4.69 


< 


T 


< 


7.12 yr 


HD 74156 c 


5.40 


< 




, < 


7.31 


5.85 


< 


T 


< 


6.38 yr 


HD 38529 c 


1.17 


< 


Xv 




0.86 


5.93 


< 


T 


< 


6.08 yr 


55 Cnc d 


19.46 


< 


Xi, 


. < 


26.33 


12.99 


< 


T 


< 


16.25 


yr 


47 UMa c 


0.66 


< 


Xv 


;< 


0.90 


5.91 


< 


T 


< 


9.22 yr 


Jupiter (2 pc) 


1.07 


< 


Xi, 


;< 


1.45 


11.11 


< 


T 


< 


11.61 


yr 


Jupiter (3 pc) 


0.87 


< 


Xp 


, < 


1.18 


10.98 


< 


T 


< 


11.66 


yr 


Jupiter (10 pc) 


0.70 


< 


Xi, 


. < 


0.95 


10.16 


< 


T 


< 


12.07 


yr 
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Table 4. Ratios of the rms errors a m for a given parameter (M p , a, T, e, and i) derived from a 

multiple-planet fit (in simulations where the full planetary system was generated around the 
parent star) to the rms errors a s derived from a single-planet solution (in simulations where only 



one planet orbited the central star). 


mi l 

. The values 


are averaged over 


il ■ l - i • 

the inclination 


angle 


Planet 


a m (M p )/a s (M p ) a m 


(a)/a s (a) , 


? m (T)/a s (T) a m 


(e)/a s (e) a m 


(*)/*-(*) 


v And b 


o /in 


2.32 


1.56 


0.98 


1.36 


v And c 


1 AO 


1.35 


1.59 


1.40 


1.35 


v And d 


1.O0 


1.37 


1.54 


1.52 


1.35 


oo One D 


i.oy 


1.39 


1.49 


1.35 


1.23 


55 Cnc c 


1 A Q 
1.40 


1.48 


1.40 


1.14 


1.18 


55 Cnc d 


1 nn 
l.UU 


0.76 


0.83 


1.53 


1.71 


47 UMa b 


15.4 


18.0 


9.64 


1.54 


2.31 


47 UMa c 


1.04 


1.09 


1.12 


3.69 


2.71 


Gliese 876 b 


1.21 


1 91 


1 91 


1 3Q 


1 19 


Gliese 876 c 


1.22 


1.23 


1.55 


1.14 


1.12 


HD 12661 b 


1.15 


1.15 


1.24 


1.15 


1.12 


HD 12661 c 


1.19 


1.21 


1.37 


1.33 


1.10 


HD 160691 b 


1.54 


1.54 


1.86 


1.38 


1.32 


HD 160691 c 


2.09 


2.07 


2.00 


1.78 


1.66 


HD 168443 b 


1.13 


1.15 


1.21 


1.13 


1.09 


HD 168443 c 


0.46 


0.47 


0.86 


0.97 


1.09 


HD 37124 b 


1.18 


1.17 


1.27 


1.16 


1.11 


HD 37124 c 


0.28 


0.48 


1.00 


0.84 


0.86 


HD 38529 b 


1.19 


1.19 


1.19 


1.04 


1.07 


HD 38529 c 


0.86 


0.85 


0.78 


0.55 


0.87 


HD 74156 b 


1.08 


1.08 


1.29 


1.06 


1.06 


HD 74156 c 


0.59 


0.71 


0.94 


0.68 


0.66 


HD 82943 b 


2.87 


2.08 


1.69 


3.23 


1.77 


HD 82943 c 


2.50 


1.81 


1.95 


2.05 


1.70 



-47- 



Table 5. Results from a fit to the fractional errors on orbital elements and planet mass, 
assuming a three-dimensional power-law dependence on the three parameters a/ad, T, and e of 
the form: 770 x (a/a d ) Vl x T m x e m . The results are averaged over the inclination angle 



Parameter Error (%) 


Best-Fi 


t Power-Law 






34.96 x (a/a d ) 


-0.72 x ^-0.42 x 


g 0.35 


0"a 


34.08 x (a/a d ) 


-0.73 x ji-0.43 x 


e 0.36 


<7t 


2.50 x (a/a d y 


-1.19 x T 2.65 x e - 


-0.12 


0"e 


23.39 x (a/a d ) 


-0.87 x T 0.01 x e 


-1.06 


(?\ 


24.28 x (a/a d ) 


-0.59 x ^-0.13 x 


e 0.08 




67.85 x (a/a d ) 


-0.43 x ^-0.06 x 


e 0.43 
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